# ===============================================================#
#                     Replication files for:                     #
#.  "Attitudinal and Behavioral Legacies of Wartime Violence:    #
#                      A Meta-Analysis"                          #
#                        Joan Barceló                            #
#               American Political Science Review                #
#               Last update: September 3, 2025                   #
# ===============================================================#

#################################
# Appendix H.1 / Table H1: Publication-bias adjusted effects
#################################

## ------------------------ Load Packages ------------------------

library(dplyr)
library(tidyr)
library(meta)
library(knitr)
library(kableExtra)
library(metasens)

# ---------- Settings ----------
min_n <- 10   # minimum number of UNIQUE studies required to compute Rücker’s limit correction

# ---------- Paths and file maps ----------
BASE <- "~/Datasets/"

read_meta <- function(path) read.csv(file.path(BASE, path), sep = ",", stringsAsFactors = FALSE)

# Outcome -> filename maps (Full)
files_full <- c(
  "Social groups participation" = "meta.groups.csv",
  "Political participation"     = "meta.political_participation.csv",
  "Voting"                      = "meta.voting.csv",
  "Generalized trust"           = "meta.generalized_trust.csv",
  "Interest / knowledge"        = "meta.political_interest.csv",
  "Normative prosociality"      = "meta.normative.csv",
  "Community leadership"        = "meta.leadership.csv",
  "Altruism"                    = "meta.altruism.csv",
  
  "Attitudes toward wartime enemies" = "meta.warenemies.csv",
  "Threat perceptions"               = "meta.threat.csv",
  "Intergroup attitudes"             = "meta.intergroup.csv",
  
  "Ingroup trust"        = "meta.ingroup_trust.csv",
  "Group voting"         = "meta.groupvoting.csv",
  "Group identification" = "meta.groupid.csv",
  
  "Political intolerance"        = "meta.political_intolerance.csv",
  "Authoritarian attitudes"      = "meta.authoritarian.csv",
  "Hawkish security preferences" = "meta.hawkish.csv",
  "Support for punitive justice" = "meta.punitive.csv",
  "Extreme ideology"             = "meta.xtrideology.csv",
  "Social intolerance"           = "meta.social_intolerance.csv",
  "Antagonism toward peace"      = "meta.antipeace.csv",
  "Institutional distrust"       = "meta.instmistrust.csv"
)

# Derive Bivariate and Quasi-Experimental file maps
files_biv <- sub("\\.csv$", "_biv.csv", files_full)
files_qe  <- sub("\\.csv$", "_random.csv", files_full)

# Exclude QE outcomes that cannot run
files_qe  <- files_qe[!names(files_qe) %in% c("Normative prosociality", "Political intolerance")]

# Master row order for final table
row_order <- names(files_full)

## ------------------------ Load data ------------------------

datasets_full <- list()
for (nm in names(files_full)) {
  datasets_full[[nm]] <- read_meta(files_full[[nm]])
}

datasets_biv <- list()
for (nm in names(files_biv)) {
  datasets_biv[[nm]] <- read_meta(files_biv[[nm]])
}

datasets_qe <- list()
for (nm in names(files_qe)) {
  datasets_qe[[nm]] <- read_meta(files_qe[[nm]])
}

# ---------- Stars ----------
stars <- function(p) {
  ifelse(is.na(p), "",
         ifelse(p < 0.001, "***",
                ifelse(p < 0.01, "**",
                       ifelse(p < 0.05, "*", ""))))
}

# ---------- FULL models ----------
res_full <- data.frame(Outcome = names(datasets_full),
                       est = NA_real_, se = NA_real_, p = NA_real_,
                       stringsAsFactors = FALSE)

for (i in seq_len(nrow(res_full))) {
  outcome <- res_full$Outcome[i]
  df <- datasets_full[[outcome]]
  nm <- names(df)
  
  if (all(c("yi","vi") %in% nm)) {
    yi <- suppressWarnings(as.numeric(df$yi))
    vi <- suppressWarnings(as.numeric(df$vi))
  } else {
    y_col  <- intersect(c("coef","yi","TE","te","estimate","est","effect","effect_size","g","d","b"), nm)
    se_col <- intersect(c("se","seTE","stderr","std.error","std_error"), nm)
    vi_col <- intersect(c("vi","variance","var"), nm)
    if (length(y_col) && length(vi_col)) {
      yi <- suppressWarnings(as.numeric(df[[y_col[1]]]))
      vi <- suppressWarnings(as.numeric(df[[vi_col[1]]]))
    } else if (length(y_col) && length(se_col)) {
      yi <- suppressWarnings(as.numeric(df[[y_col[1]]]))
      se <- suppressWarnings(as.numeric(df[[se_col[1]]]))
      vi <- se^2
    } else {
      yi <- NA_real_; vi <- NA_real_
    }
  }
  
  ok <- is.finite(yi) & is.finite(vi) & vi > 0
  yi <- yi[ok]; vi <- vi[ok]
  
  if ("authoryear" %in% nm) {
    k_studies <- length(unique(as.character(df$authoryear[ok])))
  } else {
    k_studies <- length(yi) 
  }
  
  if (k_studies >= min_n) {
    m <- tryCatch(metagen(TE = yi, seTE = sqrt(vi), sm = "GEN"),
                  error = function(e) NULL)
    if (!is.null(m)) {
      lmobj <- tryCatch(metasens::limitmeta(m), error = function(e) NULL)
      if (!is.null(lmobj)) {
        res_full$est[i] <- as.numeric(lmobj$TE.adjust)
        res_full$se[i]  <- as.numeric(lmobj$seTE.adjust)
        res_full$p[i]   <- as.numeric(lmobj$pval.adjust)
      }
    }
  }
}

res_full$Full <- ifelse(is.na(res_full$est),
                        "-",
                        sprintf("%.2f%s (%.2f)", res_full$est, stars(res_full$p), res_full$se))
res_full <- res_full[, c("Outcome","Full")]
print(res_full)  # Column Full

# ---------- BIVARIATE models ----------
res_biv <- data.frame(Outcome = names(datasets_biv),
                      est = NA_real_, se = NA_real_, p = NA_real_,
                      stringsAsFactors = FALSE)

for (i in seq_len(nrow(res_biv))) {
  outcome <- res_biv$Outcome[i]
  df <- datasets_biv[[outcome]]
  nm <- names(df)
  
  if (all(c("yi","vi") %in% nm)) {
    yi <- suppressWarnings(as.numeric(df$yi))
    vi <- suppressWarnings(as.numeric(df$vi))
  } else {
    y_col  <- intersect(c("coef","yi","TE","te","estimate","est","effect","effect_size","g","d","b"), nm)
    se_col <- intersect(c("se","seTE","stderr","std.error","std_error"), nm)
    vi_col <- intersect(c("vi","variance","var"), nm)
    if (length(y_col) && length(vi_col)) {
      yi <- suppressWarnings(as.numeric(df[[y_col[1]]]))
      vi <- suppressWarnings(as.numeric(df[[vi_col[1]]]))
    } else if (length(y_col) && length(se_col)) {
      yi <- suppressWarnings(as.numeric(df[[y_col[1]]]))
      se <- suppressWarnings(as.numeric(df[[se_col[1]]]))
      vi <- se^2
    } else {
      yi <- NA_real_; vi <- NA_real_
    }
  }
  
  ok <- is.finite(yi) & is.finite(vi) & vi > 0
  yi <- yi[ok]; vi <- vi[ok]
  
  if ("authoryear" %in% nm) {
    k_studies <- length(unique(as.character(df$authoryear[ok])))
  } else {
    k_studies <- length(yi)
  }
  
  if (k_studies >= min_n) {
    m <- tryCatch(metagen(TE = yi, seTE = sqrt(vi), sm = "GEN"),
                  error = function(e) NULL)
    if (!is.null(m)) {
      lmobj <- tryCatch(metasens::limitmeta(m), error = function(e) NULL)
      if (!is.null(lmobj)) {
        res_biv$est[i] <- as.numeric(lmobj$TE.adjust)
        res_biv$se[i]  <- as.numeric(lmobj$seTE.adjust)
        res_biv$p[i]   <- as.numeric(lmobj$pval.adjust)
      }
    }
  }
}

res_biv$Bivariate <- ifelse(is.na(res_biv$est),
                            "-",
                            sprintf("%.2f%s (%.2f)", res_biv$est, stars(res_biv$p), res_biv$se))
res_biv <- res_biv[, c("Outcome","Bivariate")]
print(res_biv)  # Column Bivariate

# ---------- QUASI-EXPERIMENTAL models ----------
res_qe <- data.frame(Outcome = names(datasets_qe),
                     est = NA_real_, se = NA_real_, p = NA_real_,
                     stringsAsFactors = FALSE)

for (i in seq_len(nrow(res_qe))) {
  outcome <- res_qe$Outcome[i]
  df <- datasets_qe[[outcome]]
  nm <- names(df)
  
  if (all(c("yi","vi") %in% nm)) {
    yi <- suppressWarnings(as.numeric(df$yi))
    vi <- suppressWarnings(as.numeric(df$vi))
  } else {
    y_col  <- intersect(c("coef","yi","TE","te","estimate","est","effect","effect_size","g","d","b"), nm)
    se_col <- intersect(c("se","seTE","stderr","std.error","std_error"), nm)
    vi_col <- intersect(c("vi","variance","var"), nm)
    if (length(y_col) && length(vi_col)) {
      yi <- suppressWarnings(as.numeric(df[[y_col[1]]]))
      vi <- suppressWarnings(as.numeric(df[[vi_col[1]]]))
    } else if (length(y_col) && length(se_col)) {
      yi <- suppressWarnings(as.numeric(df[[y_col[1]]]))
      se <- suppressWarnings(as.numeric(df[[se_col[1]]]))
      vi <- se^2
    } else {
      yi <- NA_real_; vi <- NA_real_
    }
  }
  
  ok <- is.finite(yi) & is.finite(vi) & vi > 0
  yi <- yi[ok]; vi <- vi[ok]
  
  if ("authoryear" %in% nm) {
    k_studies <- length(unique(as.character(df$authoryear[ok])))
  } else {
    k_studies <- length(yi)
  }
  
  if (k_studies >= min_n) {
    m <- tryCatch(metagen(TE = yi, seTE = sqrt(vi), sm = "GEN"),
                  error = function(e) NULL)
    if (!is.null(m)) {
      lmobj <- tryCatch(metasens::limitmeta(m), error = function(e) NULL)
      if (!is.null(lmobj)) {
        res_qe$est[i] <- as.numeric(lmobj$TE.adjust)
        res_qe$se[i]  <- as.numeric(lmobj$seTE.adjust)
        res_qe$p[i]   <- as.numeric(lmobj$pval.adjust)
      }
    }
  }
}

res_qe$`Quasi-Experimental` <- ifelse(is.na(res_qe$est),
                                      "-",
                                      sprintf("%.2f%s (%.2f)", res_qe$est, stars(res_qe$p), res_qe$se))
res_qe <- res_qe[, c("Outcome","Quasi-Experimental")]
print(res_qe)  # Column Quasi-experimental
